Two-dimensional infrared-Raman spectroscopy as a probe of water’s tetrahedrality

Two-dimensional spectroscopic techniques combining terahertz (THz), infrared (IR), and visible pulses offer a wealth of information about coupling among vibrational modes in molecular liquids, thus providing a promising probe of their local structure. However, the capabilities of these spectroscopies are still largely unexplored due to experimental limitations and inherently weak nonlinear signals. Here, through a combination of equilibrium-nonequilibrium molecular dynamics (MD) and a tailored spectrum decomposition scheme, we identify a relationship between the tetrahedral order of liquid water and its two-dimensional IR-IR-Raman (IIR) spectrum. The structure-spectrum relationship can explain the temperature dependence of the spectral features corresponding to the anharmonic coupling between low-frequency intermolecular and high-frequency intramolecular vibrational modes of water. In light of these results, we propose new experiments and discuss the implications for the study of tetrahedrality of liquid water.

Two-dimensional infrared-Raman spectroscopy as a probe of water's tetrahedrality Tomislav Begušić 1 & Geoffrey A. Blake 1,2 Two-dimensional spectroscopic techniques combining terahertz (THz), infrared (IR), and visible pulses offer a wealth of information about coupling among vibrational modes in molecular liquids, thus providing a promising probe of their local structure. However, the capabilities of these spectroscopies are still largely unexplored due to experimental limitations and inherently weak nonlinear signals. Here, through a combination of equilibrium-nonequilibrium molecular dynamics (MD) and a tailored spectrum decomposition scheme, we identify a relationship between the tetrahedral order of liquid water and its two-dimensional IR-IR-Raman (IIR) spectrum. The structure-spectrum relationship can explain the temperature dependence of the spectral features corresponding to the anharmonic coupling between low-frequency intermolecular and high-frequency intramolecular vibrational modes of water. In light of these results, we propose new experiments and discuss the implications for the study of tetrahedrality of liquid water.
Understanding the dynamical structure of liquid water is paramount to a number of chemical and biological processes. In particular, the tetrahedral ordering of molecules, stemming from the directionality of hydrogen bonds, has been proposed as the origin of water's anomalous behavior in the liquid phase 1 . However, the tetrahedral structure of water has also been contested both computationally 2 and experimentally 3 . To date, most of our microscopic, structural information about liquid water comes from molecular dynamics (MD) simulations, which depend strongly on the choice of electronic structure theory or force field parametrization. In contrast, experimental tools capable of studying the local arrangement of water molecules remain scarce 4 . For example, common techniques that probe the structure of liquids, such as X-ray and neutron scattering, typically report on highly time-averaged quantities and can be accurately reproduced with very different structural motifs 2,[5][6][7] . Vibrational spectroscopy, such as IR absorption and Raman scattering, offers complementary information about the strength of intermolecular hydrogen bonds, which depends on instantaneous, local arrangement of water molecules 8 . Different spectral regions have been probed to investigate intermolecular hydrogen-bond bending and stretching modes (up to 300 cm −1 ), frustrated rotational (librational) modes (400−1000 cm −1 ), intramolecular bending (1650 cm −1 ) and stretching (3000−3800 cm −1 ) modes, as well as combination bands between intermolecular and intramolecular modes. Reported experiments and simulations range from high-resolution spectroscopy on molecular clusters [9][10][11][12] to studies of interfacial [13][14][15] or bulk water [16][17][18][19][20][21][22][23][24][25] . Recently, a joint experimental and computational study 26,27 of the temperature dependence of the Raman spectrum of liquid water revealed a clear structure-spectrum relationship between the local tetrahedral order parameter 28 , a measure of structuring of liquid water, and frequency shifts and intensities of spectral peaks. In fact, the temperature dependence of the spectrum was, in this way, fully explained as a change in the thermal distribution of the tetrahedral order parameter. Yet, conventional, steady-state spectroscopy of liquids typically produces broad, unresolved spectral features and misses dynamical information.
Time-resolved and two-dimensional vibrational spectroscopies have emerged in the past years as probes sensitive to local water coordination. Even here, however, the long-standing tool of twodimensional infrared spectroscopy (2DIR) 4,29,30 , which measures the interaction of light with intramolecular, high-frequency modes, provides only an indirect probe of intermolecular dynamics. 2DIR has been successfully applied to analyze coupling among high-frequency modes in ice 31 and in proteins 32 , but also to understand the vibrational relaxation mechanisms in liquids 33 . For example, several 2DIR studies on liquid water and ice speculated that mechanical anharmonic coupling to intermolecular modes contributes to the relaxation of the OH stretch [34][35][36] . In addition, non-Condon effects, which are related to electrical anharmonic coupling between high-and low-frequency modes, were shown to affect the 2DIR spectra of OH stretching mode in liquid water 37 . To target the low-frequency modes directly, a number of hybrid spectroscopic techniques have been proposed, involving different sequences of THz, IR, and visible pulses, such as the THz-THz-Raman 38-41 , THz-Raman-THz, Raman-THz-THz [42][43][44][45] , and THz-IR-Raman (also called THz-IR-visible 46,47 or TIRV). Their development was enabled by the recent advances in the generation of strong THz pulses that are needed to induce a nonlinear light-matter interaction 48 . The THz-THz-Raman and TIRV methods are related to other two-dimensional IR-Raman techniques, namely the two-dimensional IR doubly vibrationally enhanced and IR-IR-visible sum-frequency generation spectroscopies [49][50][51][52][53][54] . For example, TIRV experiments have revealed unambiguous spectral signatures of coupling between the intramolecular O-H stretch and intermolecular hydrogen-bond bending and stretching modes 46 . Theoretical simulations by Ito and Tanimura 55 predicted such spectral features and assigned them to both mechanical and electrical anharmonic coupling between the said vibrational modes. Similarly, THz-THz-Raman spectroscopy recently revealed signatures of anharmonic coupling between phonons of ionic solid LiNbO 3 56 . Finally, Raman-THz-THz and THz-Raman-THz spectroscopies have been used to study the inhomogeneity of liquid water and aqueous solutions 57,58 , as well as the coupling among intermolecular and intramolecular modes in liquid and solid bromoform 59,60 . Even so, due to limited availability of efficient THz emitter materials, not all frequencies have been covered by the reported techniques. Specifically, most two-dimensional hybrid THz-Raman spectroscopies of liquid water targeted hydrogen-bond bending and stretching modes, i.e., frequencies up to 400 cm −1 , leaving the water librational dynamics largely unexplored.
Here, we aim to provide new insights into the capabilities of twodimensional hybrid IR-Raman vibrational spectroscopies. To this end, we study the temperature dependence of the two-dimensional IR-IR-Raman (IIR) spectrum, which is given by the double Fourier (or sine) transform of an appropriate third-order, two-time response function. Since the computational model involves all vibrational modes of the system, the response function covers a broad range of frequencies, which, in practice, can be mapped out only through separate experiments. To date, only the TIRV frequency region has been experimentally measured, although its temperature dependence was not studied. Second, following ref. 26, we then establish a structure-spectrum relationship by separating the spectral contributions from molecules exhibiting low or high tetrahedral coordination. Third, to justify the molecular dynamics (MD) results at low temperatures, we analyze whether nuclear quantum effects are discernible in the twodimensional IIR spectrum.

Theoretical model
We simulated the IIR response function 40,55 Rðt 1 ,t 2 Þ = À 1 _ 2 Trf½½Πðt 1 + t 2 Þ,μðt 1 Þ,μð0Þρg ð1Þ using the equilibrium-nonequilibrium MD approach, in which the quantum-mechanical trace is replaced by a classical average 61-64 whereρ is the thermal density operator,μ = μðqÞ is the dipole moment operator, andΠ = ΠðqÞ is the polarizability. q t denotes the position vector of a classical trajectory initiated at (q 0 , q 0 ), whereas q ±,t corresponds to the initial conditions ðq ± ,0 = q 0 ,p ± ,0 = p 0 ± εμ 0 ðq 0 Þ=2Þ after an instantaneous interaction with the electric field. ε is a free parameter in the calculations and corresponds to the magnitude of the external electric field integrated over the short interaction time. For sufficiently small ε, the thermal average of Eq. (2) is linear in ε, meaning that the time-dependent response function R MD (t 1 , t 2 ) is, as expected, independent of its exact value 63,65 . The two-dimensional spectra are computed through a double sine transform 42,55 Rðω 1 ,ω 2 Þ = To model water, we used a flexible, point-charge qTIP4P/F force field 66 , which has been well studied for spectroscopic simulations 25 and benchmarked against a number of experimental thermodynamic properties of water, including the radial distribution functions, dielectric constant, density, and melting point. As a point-charge force field, the model is computationally efficient, which is needed for statistically averaging two-time response functions. In addition, because it is flexible, we could study all vibrational degrees of freedom, including high-frequency intramolecular modes. To allow for nonlinear dependence of the dipole moments and polarizabilities on nuclear coordinates, we employed the truncated dipole-induceddipole (DID) model. Following Hamm 44 , each water molecule was amended with permanent anisotropic polarizability, which was used for the evaluation of the induced contributions to the dipoles and polarizabilities. In contrast to the rigid water simulations of ref. 44, here the permanent polarizability was an explicit function of intramolecular degrees of freedom, according to ref. 67. The introduced induced dipole and polarizability effects were not used in the evaluation of the forces, which were computed according to the original, non-polarizable qTIP4P/F model. We note that more advanced models for the potential energy, dipoles, and polarizabilities of liquid water exist and have been used in the simulation of one-and twodimensional spectroscopies 19,22,24,26,44,68 . However, while some of them were not readily available, others considered only rigid water molecules or were fitted to experiments using classical MD simulations, which would prevent us from accurately exploring nuclear quantum effects. In Fig. 1, we show that the employed force field combined with our induced dipole and polarizability models can reproduce the main features of the experimental IR absorption and anisotropic Raman spectra of liquid water. Additional details about the model and MD simulations can be found in the Methods section.
Most of the results presented below rely on the validity of the classical MD approach, which neglects the quantum-mechanical properties of atomic nuclei. However, due to the presence of light hydrogen atoms, their effect might not be negligible 66 . Although such nuclear quantum effects on one-dimensional IR and Raman spectra have been well studied using approximate but reliable classical-like methods, no tools similar to the equilibrium-nonequilibrium MD have been available to study nuclear quantum effects on two-dimensional IR-Raman spectra [69][70][71] . Recently, we have developed a new ringpolymer MD (RPMD) approach 65 , which can simulate, at least approximately, the nuclear quantum effects on the two-dimensional IIR spectra. Briefly, the RPMD method 72 replaces the original quantummechanical problem by an extended classical system consisting of N replicas (beads) of the original system connected by harmonic springs. For a given potential energy surface, the extended classical system in the limit of N → ∞ reproduces the exact quantum-mechanical thermal distribution of nuclear degrees of freedom, while if N = 1, RPMD reduces to classical MD. In our simulations, we used N = 32, which is sufficiently large for liquid water in the studied temperature range 25 . Since RPMD is known to suffer from the spurious resonance issue, where unphysical peaks due to artificial harmonic springs appear in the spectra, we employed its thermostatted version (TRPMD) 73,74 . IR absorption and anisotropic Raman spectra simulated with TRPMD are presented in Fig. 1. Figure 2 shows the full two-dimensional IIR spectrum of liquid water simulated at 300 K. Apart from the spectral features along the ω 1 = ω 2 diagonal, which appear due to mechanical or electrical anharmonicity of individual vibrational modes, the spectrum contains off-diagonal peaks that correspond to coupling among different vibrations. Overall, the spectrum qualitatively agrees with the simulation of ref. 55, with the main difference in the relative intensities of different spectral regions, which can depend strongly on the details of the potential energy, dipole, and polarizability surfaces. In the following, we focus on the two frequency regions indicated by the grey rectangles. One region targets the coupling between intermolecular modes (0 cm −1 < ω 1 < 1200 cm −1 ) and the intramolecular O-H stretch mode (2900 cm −1 < ω 2 < 4100 cm −1 ), whereas the other covers all low-frequency, intermolecular modes, and the intramolecular bend mode.  Fig. 3) imply that the complex spectral lineshape of the former, which we will refer to as the TIRV region 46 , are a product of interplay between mechanical and electrical anharmonicity. Namely, mechanical anharmonicity leads to a shape comprising a positive and a negative lobe above and below the central ω 2 frequency, whereas electrical anharmonicity produces an approximately symmetric feature. To verify this interpretation, we would ideally construct approximate dipole and polarizability models that depend linearly on atomic coordinates, which would allow us to study the mechanical anharmonic coupling pathways directly. This can be easily done for the dipole moments by neglecting the induced part because the permanent molecular dipole is a linear function of coordinates (see Methods). Unfortunately, the same cannot be achieved for polarizability because both permanent and induced parts are nonlinear in atomic coordinates. For this reason, we consider an alternative, Raman-IR-IR (RII) pulse sequence, in which the polarizability is responsible for the first interaction with the external electric field. Assuming weak electrical and mechanical anharmonic coupling, it can be shown that the nonlinearity in the first interaction does not contribute to the RII spectrum (Supplementary Discussion 2). Therefore, the RII spectrum with permanent (i.e., linear) dipole moments (shown in Fig. 3a) consists only of mechanical anharmonic coupling pathways. Because the Raman response in the THz frequency range is weak, the corresponding RII spectrum exhibits roughly an order of magnitude lower intensity than the TIRV spectrum. For comparison with the full TIRV spectrum, the simulated RII spectrum must be appropriately scaled (Fig. 3b) along the frequency axes (Supplementary Eq. (9)). The result agrees with the proposed interpretation that the nodal shape of the spectrum results from mechanical anharmonic coupling pathways.

Two-dimensional IIR spectrum of liquid water
In the low-frequency region of interest, we observe a strong peak at about (600 cm −1 , 1650 cm −1 ) due to anharmonic coupling between librations (400-1000 cm −1 ) and the intramolecular bending mode 75 . The same coupling mechanism is responsible for the combination transition at around 2150 cm −1 , also known as the "association" band 9,20 , in the one-dimensional spectra. Interestingly, this combination band is not captured in our model (Fig. 1), even though the corresponding two-dimensional spectral feature clearly appears in the IIR spectrum. This discrepancy can be explained by the fact that the fundamentals of the one-dimensional spectra follow harmonic selection rules, while the peaks in the two-dimensional spectra appear solely due to the anharmonic excitation pathways. Therefore, the combination bands in the one-dimensional spectra can be orders of magnitude weaker than the fundamentals and still exhibit strong off-diagonal peaks in the two-dimensional spectra. Indeed, we see that the diagonal Fig. 2 | Two-dimensional IIR spectrum of liquid water simulated with classical MD at 300 K. Grey squares indicate the regions of the spectrum that we studied in more detail, namely the region of TIRV spectroscopy and the low-frequency region that covers intermolecular modes and the intramolecular bending mode. peak at around (1650 cm −1 , 1650 cm −1 ) is much lower in intensity than the off-diagonal libration-bending peak, implying that the anharmonicity within the bending mode is weaker than its anharmonic coupling to the intermolecular librations. Let us note that other force fields and dipole/polarizability models also heavily underestimate or fail to reproduce the combination band at 2150 cm −119,21,55,76 , unlike the ab initio approaches, which appear to systematically reproduce it 26,77,78 . Since the force field we used reproduces the structural and dynamical properties of liquid water rather accurately 66 , we tentatively assign the absence of this combination band in the simulated spectra (Fig. 1) to the limitations of our dipole and polarizability models.

Temperature dependence of the IIR spectrum
We now turn to the temperature dependence of the two IIR spectral regions (Fig. 4). The TIRV part of the spectrum experiences several changes as the temperature is increased from 280 K to 360 K. At higher temperatures, the peaks are blue-shifted along ω 2 , which agrees with the observed trends in the experimental Raman spectra and simulated vibrational density of states 26 . In addition, the shape of the spectral peaks changes drastically. The most prominent change in the spectrum is the disappearance of the negative peak around (250 cm −1 , 3600 cm −1 ). This indicates that the electrical anharmonicity contribution to the spectrum, marked by a strongly symmetric lineshape, increases compared to that of the mechanical anharmonicity. Unlike the frequency shift, this feature cannot be accessed through steady-state spectroscopy. Finally, we note that a peak around (400 cm −1 , 3700 cm −1 ) appears at elevated temperature, while the other peak at about (700 cm −1 , 3400 cm −1 ) almost disappears. Interestingly, both of these coupling terms have been discussed as possible origins of the libration-stretch combination band appearing in the anisotropic Raman spectrum of liquid water at 4100 cm −126,79 . The same combination band has also been shown to play an important role in second-order vibrational sum-frequency spectroscopy of interfacial water 15 .
The low-frequency region also exhibits strong temperature dependence. For example, the results (see also difference spectra in Supplementary Fig. 5) clearly indicate a red shift of the libration-bend peak along ω 1 with increasing temperature. This is a straightforward consequence of the equivalent frequency shift found in the linear IR absorption spectrum ( Supplementary Fig. 6), which also agrees with the experimentally observed trend 17 . Furthermore, the twodimensional IIR spectrum at 280 K contains peaks close to diagonal, around ω 1 = ω 2 = 250 cm −1 , due to anharmonicity of the hydrogen-bond stretching modes. These seem to progressively disappear at elevated temperatures. Importantly, this change in intensity could not be simply predicted from the one-dimensional spectra, which show little change in the intensity of the 250 cm −1 band. We note, however, that the features in this congested spectral region could be affected by the short times available from our simulations and by artificial broadening. Longer simulations are possible with rigid water models, which have been used to study explicitly the intermolecular modes and corresponding two-dimensional THz-Raman signals in the time domain 44,62 .
The TIRV spectral region has been studied experimentally 46,47 . However, the experiments could only measure an absolute value Fourier transform spectrum, which is quite different from the real sine transform discussed here and in ref. 55. Specifically, the experimental spectrum reported in Fig. 8     response function bỹ Sðω 1 ,ω 2 Þ = j½Sðω 1 ,ω 2 Þ + Sðω 2 À ω 1 ,ω 2 ÞE THz ðω 1 ÞE IR ðω 2 À ω 1 Þj, ð4Þ where E THz (ω) and E IR (ω) are the THz and IR electric fields obtained as square roots of the pulse intensity spectra presented in Fig. 2b, c of ref. 46. The main differences between our simulation and the experiment of ref. 47 arise due to the limitations of our model, namely the narrow and blue-shifted OH stretch (Fig. 1). Furthermore, the experimental spectrum exhibits additional spectral features that cannot be explained with Eq. (4) and the pulse shapes we used because they fall outside of the bandwidth of the instrument response function. Unlike the sine-transform spectra shown in Fig. 4, Fouriertransform spectra convolved with external electric fields (Fig. 5) only become less intense with increasing temperature but show no interesting spectral change. Therefore, an accurate determination of the full response function will be needed to experimentally measure the spectral features appearing in the sine-transform TIRV spectra. Fortunately, a scheme that could achieve this goal has been recently applied to the TIRV spectrum of liquid dimethyl sulfoxide 80 , demonstrating that similar experiments on water are within reach.

Spectral signatures of water's tetrahedrality in the IIR spectrum
To understand the temperature-dependent spectral features, we follow the work of Morawietz et al. 26 , where the temperature effects were analyzed in relation to the local structuring around individual water molecules. This can be quantified by the local tetrahedral order parameter Q 81-83 , which measures how the arrangement of the four neighboring water molecules deviates from the ideal tetrahedral arrangement around the central one. By convention, Q = 1 for a perfect tetrahedral arrangement, while Q = 0 for an ideal gas. The distribution of Q at thermal equilibrium exhibits a strong temperature dependence and a clear isosbestic point at Q ≈ 0.67 (Fig. 6a). In the studied temperature range, the bimodal distribution can be decomposed into two temperature-independent components whose populations change with temperature 26 . However, as shown by Geissler 84 , the appearance of an isosbestic point does not necessarily imply that water is a heterogeneous mixture. More recently, the increased population of low-Q water molecules at higher temperatures has been assigned to the appearance of neighboring molecules in the interstitial position between the first and second solvation shells 83 .
To directly correlate the local tetrahedral order parameter with the spectral features observed in the two-dimensional IIR spectra, we decomposed the spectrum into contributions of individual molecules (Supplementary Discussion 3). Then, we simulated the spectra (Fig. 6b-e) originating predominantly from the molecules with either high (Q > 0.72) or low (Q < 0.62) local tetrahedral order parameters. The two IIR spectra align remarkably well with the observed temperature-dependent changes. Specifically, the high-order TIRV spectrum (Fig. 6c) exhibits a clear mechanical anharmonic coupling feature with a positive and a negative lobe at ω 1 ≈ 250 cm −1 , whereas the corresponding low-order spectrum contains a strong symmetric peak, indicating electrical anharmonic coupling between the hydrogen-bond modes and O-H stretch. In addition, the libration-stretching peak at (400 cm −1 , 3700 cm −1 ) appears almost exclusively in molecules with low tetrahedral order. Similarly, the libration-bending peak in the lowfrequency part the spectrum (Fig. 6d, e) appears at lower ω 1 frequencies for low Q. Overall, the temperature dependence of the IIR spectra can be almost exclusively assigned to the changes in the distribution of the local tetrahedral order parameter and the effect of the local structure on the spectral features. Therefore, TIRV spectroscopy expanded into the water libration frequency range, i.e., with ω 1 covering up to ≈ 800 cm −1 , could provide an alternative probe of the local molecular structure in liquid water and aqueous solutions. Importantly, the spectral changes in IIR are more drastic, and therefore more sensitive to local ordering, than the frequency shifts observed in conventional, one-dimensional IR and Raman spectroscopies.
This result is consistent with other experimental observations and theoretical models. Namely, it is known that as the tetrahedral order parameter increases, the hydrogen-bond stretching mode frequency increases and the OH stretching frequency decreases 26 , which can be interpreted in terms of a growing mechanical anharmonic coupling strength between the two modes. A more detailed analysis is possible based on the work of Auer and Skinner 85 , who studied the dependence of the intramolecular stretching frequency and the corresponding dipole moment derivative on the electric field E generated by the surrounding water molecules at the hydrogen atom of the central molecule and projected along the OH bond. They found that the frequency can be fit to a quadratic function of this electric field, while the dipole moment derivative is approximately linear in E: where all parameters ω ðαÞ OH and μ 0ðαÞ OH are positive and defined in Table I. of ref. 85. More recent models 86,87 included a very weak quadratic term for the dipole moment as well, which can be neglected within the typical range of values of the electric field. For us, the relevant anharmonic coupling quantities are the derivatives of ω OH and μ 0 OH Fig. 6 | Effect of tetrahedral order on the IIR spectra. a Distribution of the local tetrahedral order parameter of liquid water as a function of temperature. The isosbestic value is indicated by the solid, grey line. Two dashed, grey lines specify the region excluded from subsequent two-dimensional IIR spectra simulations (0.62 < Q < 0.72). Two-dimensional IIR spectra of water molecules with high (Q > 0.72, c, e) or low (Q < 0.62, b, d) tetrahedral order parameter, simulated with the classical equilibrium-nonequilibrium MD approach at 320 K.
with respect to an intermolecular hydrogen-bond mode q HB , ∂ω OH =∂q HB = À ðω ð1Þ OH + 2ω ð2Þ OH EÞ∂E=∂q HB and ∂μ 0 OH =∂q HB = μ 0ð1Þ OH ∂E=∂q HB . From these equations, we see that the mechanical anharmonic coupling, determined by |∂ω OH /∂q HB |, grows with the electric field, whereas the electrical anharmonic coupling, proportional to |∂μ 0 OH =∂q HB |, has no explicit dependence on E. In ref. 83, it has been shown that E correlates positively with the tetrahedral order parameter Q, which further validates the interpretation of our two-dimensional spectroscopy simulations, i.e., stronger mechanical anharmonic coupling in high-Q molecules.

Nuclear quantum effects on the two-dimensional IIR spectrum
Finally, we report the first TRPMD simulation of the two-dimensional IIR spectrum of liquid water and compare it with the MD simulation in Fig. 7. Some nuclear quantum effects can be easily explained in terms of the differences between one-dimensional spectra ( Supplementary  Fig. 7). Namely, the red-shifted O-H stretch band in the TRPMD steadystate spectra reflects in a red shift along ω 2 of the spectral features in the TIRV spectrum (Fig. 7a, b). Furthermore, the libration-stretch peak shifts from ω 1 ≈ 700 cm −1 in the MD spectrum to ω 1 ≈ 800 cm −1 in TRPMD, which aligns with the differences between MD and TRPMD IR absorption spectra. In the low-frequency part of the spectrum (Fig. 7c,  d), the libration-bend peak is analogously shifted along ω 1 . However, some features cannot be understood from comparison with the onedimensional spectra, such as the absence of the strong negative feature at (ω 1 ≈ 750 cm −1 , ω 2 ≈ 50 cm −1 ) in TRPMD or the change in the spectral lineshapes in the TIRV region at ω 1 > 400 cm −1 . Nevertheless, these differences are not sufficiently large to alter the above conclusions based on classical MD simulations.

Discussion
To conclude, we have presented MD and TRPMD simulations of the two-dimensional IIR spectra of liquid water. By analyzing the temperature dependence of the spectrum, we have identified features that report on the degree of local tetrahedral ordering in the first coordination shell. Further computational work is needed to confirm whether IIR spectra of related systems, such as aqueous solutions, alcohols, or ice, can be broadly mapped to the local tetrahedral order parameter. Our work demonstrates that the temperature dependence of such two-dimensional spectra can be used in combination with computational modeling to gain insight into the structure of complex liquids. Furthermore, the simulations presented here imply that the contributions of mechanical and electrical anharmonic coupling change as a function of local tetrahedrality, which cannot be studied with conventional, one-dimensional spectroscopy. Specifically, the electrical anharmonicity dominates at higher temperatures, where the tetrahedral order parameter is low, whereas the mechanical coupling, due to the anharmonic terms in the potential energy surface, is pronounced at lower temperature and higher tetrahedrality. This finding was related to the electric field caused by the surrounding molecules and its positive correlation to the tetrahedral order parameter. Overall, the observed change in the mechanical anharmonic coupling could impact our understanding of the OH stretch relaxation mechanism. Although this relaxation has been largely assigned to the coupling with the intramolecular bending overtone, the intermolecular, hydrogenbond modes are also believed to play an important role [34][35][36]88 . Strong mechanical anharmonic coupling between intramolecular stretching and intermolecular hydrogen-bond stretching modes could shorten the excited-state lifetime of the OH stretch. Such correlation would also agree with the experimentally observed trend of shorter OH stretch excitation lifetime at a lower temperature.
Finally, we observe that the 400 cm −1 < ω 1 < 1000 cm −1 region of the IIR spectrum, corresponding to the librations of water molecules, carries rich information about liquid water and should be further explored experimentally. In particular, such studies could provide additional information on the combination bands appearing in IR absorption and Raman scattering spectra that have challenged physical chemists for decades.

Dipole moment and polarizability models
The induced dipole moments and polarizabilities were modeled as 44 where E ij = P a2fM,H1,H2g Q ja r ja,iO =r 3 ja,iO is the electric field produced by molecule j on the oxygen atom of molecule i, T ij = T(r jO,iO ), T(r) = (r 2 − 3r ⊗ r)/r 5 is the 3 × 3 dipole-dipole interaction tensor, r ja,iO = r ja − r iO , r ja denotes the three-dimensional position vector of atom a ∈ {M, O, H1, H2} of molecule j, Q ja is its partial qTIP4P/F charge, r = ∥r∥ denotes the norm of vector r, and ⊗ is the outer product of two vectors. M denotes the M site of the TIP4P model, whose position is related to the molecular geometry of molecule i by r iM = γr iO + (1 − γ) (r iH1 + r iH2 )/2, where γ = 0.73612 66 . Permanent dipole moment of molecule i was μ i = ∑ a∈{M,H1,H2} Q ia r ia /c, where c = 1.3 is a scaling factor that reduces the dipole moment of water to its gas-phase value 44,63 . The permanent polarizability of each molecule i, Π i ≡ Π i (r iO , r iH1 , r iH2 ), was computed according to ref. 67 in the molecular reference frame and rotated into the laboratory frame. Then, the total dipole moment and polarizability are given by Finally, we note that the sums over i in all expressions above are taken over the molecules in a unit cell, while the sums over j extend to infinity beyond the central unit cell. Due to the slow convergence of Coulomb interactions in the direct space, we used Ewald summation for the electric field E ij and tensor T ij 22,63 .

Computational details
The classical thermal average of Eq. (2) in the main text was evaluated by sampling from 5 independent NVT trajectories, each equilibrated for 100 ps, of a water box with 64 molecules and the cell parameter adjusted to the experimental density at any given temperature. Initial structures were taken from ref. 24. A time step of Δt = 0.25 fs, second order symplectic integrator, and a Langevin thermostat with the time constant τ = 100 fs were used throughout. From 2560 initial samples collected in this way, we launched 25 ps NVE equilibrium trajectories, along which we collected the dipole moments and polarizabilities every 4 steps (1 fs) and the derivatives of the dipole moment every 1000 steps (250 fs). Finally, starting from positions at which the dipole derivatives were evaluated, "nonequilibrium" trajectories, i.e., with modified initial momenta p ±,0 , were propagated for 1000 steps (250 fs). This generated 2.56 × 10 5 samples for Eq. (2) of the main text, which is comparable to the numbers used in Refs. 46,62,63, and the response function was computed for 250 fs in both t 1 and t 2 . The statistical error was estimated using bootstrapping and is shown in Supplementary Fig. 1. We set ε = 2, which converts into the electric field magnitude of E 0 = ε/Δt ≈ 10 V/Å. A weaker electric field was also tested (ε = 0.2a.u., see Supplementary Fig. 2) in a simulation with the doubled number of initial conditions. RPMD simulations were performed in a similar way, using 4 independent thermostatted path-integral MD (PIMD) trajectories to produce 512 samples from which TRPMD trajectories were launched. Each TRPMD equilibrium trajectory was propagated for 50 ps (200000 steps) and nonequilibrium trajectories was launched every 1000 steps. In total, this resulted in 1.024 × 10 5 initial conditions. PIMD trajectories used a path-integral Langevin equation (PILE) thermostat for the normal modes with the centroid coupled to a global velocity rescaling thermostat with τ = 100 fs 89 , whereas the TRPMD trajectories used a generalized Langevin equation (GLE) thermostat coupled to the normal mode representation of the ring polymer 74 .
Linear absorption and anisotropic Raman spectra were computed from the equilibrium trajectories according to: I Raman ðωÞ / ÀAðωÞIm where μ is the dipole moment vector of the cell, βðtÞ = ΠðtÞÀ ITr½ΠðtÞ=3, Π(t) is the polarizability tensor, I is the 3 × 3 identity matrix, and σ t = 1000 fs. A(ω) was different for the two Raman experiments shown in Fig. 1. For the broad-range Raman spectrum shown in Fig. 1b, AðωÞ = ½1 À expðÀβ_ωÞ À127 , while for the low-frequency part shown in the inset of that figure, which was experimentally measured by optical Kerr effect 16 , A(ω) = 1. Two-dimensional IIR spectra were computed from the response function according to Eq. (3) of the main text, but with a damping f ðt 1 ,t 2 Þ = expðÀðt 1 + t 2 Þ 12 =τ 12 Þ, τ 12 = 5 × 10 28

Data availability
The data generated in this study, together with the input files and scripts for plotting the data, have been deposited in the Zenodo database under accession code DOI:10.5281/zenodo.7619094.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.